#include <math.h>

double _Integral(double a, double b, double (*f)(double), double epsilon)
{
	int i, n;
	double h, x, sum;
	double t1, t2, s1, s2, c1, c2, r1, r2;
	n = 1;
	h = b - a;
	t1 = (h * ((*f)(a) + (*f)(b))) / 2;
	while(1)
	{
		sum = 0;
		for(i = 0; i < n; i++)
		{
			x = a + (h / 2) + i * h;
			sum += (*f)(x);
		}
		t2 = (t1 + sum *h) / 2;
		s2 = (4 * t2 - t1) / 3;
		if(n == 1)
		{
			t1 = t2;
			s1 = s2;
			h /= 2;
			n *= 2;
			continue;
		}
		c2 = (16 * s2 - s1) / 15;
		if(n == 2)
		{
			t1 = t2;
			s1 = s2;
			c1 = c2;
			h /= 2;
			n *= 2;
			continue;
		}
		r2 = (64 * c2 - c1) / 63;
		if(n == 4)
		{
			t1 = t2;
			s1 = s2;
			c1 = c2;
			r1 = r2;
			h /= 2;
			n *= 2;
			continue;
		}
		if(fabs(r2 - r1) < epsilon)
			break;
		t1 = t2;
		s1 = s2;
		c1 = c2;
		r1 = r2;
		h /= 2;
		n *= 2;
	}
	return r2;
}
